Solidlike behavior and anisotropy in rigid frictionless bead assemblies 



Pierre-Emmanuel PeynearO and Jean-Noel Roux 
Universite Paris-Est, UR Navier, LMSGCf , 2 allee Kepler, 
Cite Descartes, 77420 Champs-sur-Marne, France 
(Dated: Accepted for publication in Phys. Rev. E, October 1, 2008) 

We investigate the structure and mechanical behavior of assemblies of frictionless, nearly rigid 
equal-sized beads, in the quasistatic limit, by numerical simulation. Three different loading paths 
are explored: triaxial compression, triaxial extension and simple shear. Generalizing recent result, 
we show that the material, despite rather strong finite sample size effects, is able to sustain a fi- 
nite deviator stress in the macroscopic limit, along all three paths, without dilatancy. The shape 
of the yield surface in principal stress space differs somewhat from the Mohr-Coulomb prediction, 
and is more adequately described by the Lade-Duncan or Matsuoka-Nakai criteria. We study ge- 
ometric characteristics and force networks under varying stress levels within the supported range. 
Although the scalar state variables stay equal to the values observed in systems under isotropic 
pressure, the material, once subjected to a deviator stress, possesses some fabric and force distri- 
bution anisotropies. Each kind of anisotropy can be described, in good approximation, by a single 
parameter. Within the supported stress range, along each one of the three investigated stress paths, 
among those three quantities: deviator stress to mean stress ratio, fabric anisotropy parameter, 
force anisotropy parameter, any one determines the values of the two others. The pair correlation 
function also exhibits short range anisotropy, up to a distance between bead surfaces of the order 
of 10% of the diameter. The tensor of elastic moduli is shown to possess a nearly singular, uniaxial 
structure related to stress anisotropy. Possible stress-strain relations in monotonic loading paths 
are also discussed. 

PACS numbers: 45.70.-n, 83.80.Hj, 81.40.Lm, 83.10.Rs 



I. INTRODUCTION 

The disordered packing of rigid, frictionless spheri- 
cal balls epitomizes the large class of materials made 
of athermal, amorphous assemblies of particles with ex- 
tremely short range interactions, such as granular ma- 
terials P, H, H, ff [E], concentrated suspensions [1, 0], 
or some glasses Q- Obviously a highly idealized mate- 
rial, it isperhags exclusi vely studied by numerical simu- 
lations [lIlIO, Ii1j[3- However, its main merit is 
to capture the essential role of steric exclusion and pack- 
ing geometry in the rheological properties of many dif- 
ferent materials (termed "jammed" [15] in the recent lit- 
erature). In general, "jammed" particulate systems with 
strongly repulsive interactions tend to behave like plas- 
tic solids for nearly isotropic stress states, and to flow 
like liquids once the deviatoric stress reaches some fail- 
ure threshold. The solidlike regime of granular materials 
has long been described and modeled at the continuum 
scale in the field of soil mechanics 0, [l^, [l^ ■ 

Assemblies of frictionless and cohesionless grains, in 
the rigid limit, have two remarkable properties \18^ . 
First, equilibrium configurations under specified exter- 
nally applied loads minimize potential energy, thereby 
satisfying geometric optimization criteria. In particu- 
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lar, equilibrium states under isotropic pressure realize a 
local maximum of density in configuration space, sub- 
ject to impenetrability constraints. One thus obtains, 
with assembling procedures fast enough to bypass incip- 
ient crystallisation, the so-called random close packing 
(RCP) states of sphere packings [HI, [13, [13], with solid 
fraction <i> ~ 0.64. Then, the force-carrying contact net- 
work (the backbone) is generically devoid of force inde- 
terminacy, and even isostatic with circular or spherical 
objects [lH, [H, [11] ■ Consequently, equilibrium forces are 
geometrically determined, as well as the load increments 
necessary to destabilize contact networks; and such ma- 
terials, in the solid state, tend to deform in a sequence of 
rearranging events, in which the contact structure gets 
continuously broken and repaired [H, [l^, [13] . 

In spite of those appealing properties of frictionless 
spheres (or disks in 2D) , which highlight the connections 
between geometry and mechanics and endow them with 
quite generic features, the study of those model mate- 
rials is still incomplete in the published literature. Nu- 
merical investigations have mostly focused on the geom- 
etry of RCP states [H M, Il3 |, o n the possible effects of 
confining pressure variations 2l|, [2^ and specific elastic 
properties 23^ 24] , o n the one hand; and on steady-state 
shear flows [2^ |2^ [53] on the other hand. The solid 
range, in which moderate deviator stresses are supported 
by anisotropic packings in equilibrium, has hardly been 
investigated. 

In a recent publication [l3], we checked that rigid, fric- 
tionless bead packings have a finite macroscopic friction 
coefficient /i* in simple shear, and showed them to be de- 
void of dilatancy, unlike dense frictional grain assemblies. 
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The non- vanishing value of fj,* was attributed to the pos- 
sibihty to form equiUbrium structures with anisotropic 
contact networks. Both static (yield threshold) and dy- 
namic (i.e., measured in steady shear flows) values of /i* 
were shown to agree in the limit of large samples. 

The present paper further investigates the mechanical 
properties of solidlike assemblies of frictionless beads un- 
der quasistatic loading conditions. The model material, 
initially assembled under an isotropic pressure, is sub- 
jected to different deviatoric loading paths (Section [ill, 
so that a failure criterion, or yield surface, delineating 
the stable solid range in stress space can be identified in 
the macroscopic limit (Section -thus generalizing the 
friction angle measured in simple shear. Then we study 
the geometric and micromechanical features of equilib- 
rium states throughout the supported range of stresses, 
generalizing the results obtained on isotropic packings 
to systems with various levels and various directions of 
anisotropy (Section lTV]) . Stresses are related to fabric and 
force distribution anisotropy parameters by a simple for- 
mula. We show how the anomalous elastic moduli tensor 
of frictionless, nearly rigid networks is affected by stress 
anisotropy. A macroscopic stress-strain relation appears 
to be approached along the investigated monotonic load- 
ing paths, in spite of large statistical uncertainties. The 
paper ends with a discussion (Section |V|. 

II. MODEL MATERIAL AND SIMULATED 
MECHANICAL TESTS 

A. Constituents and microscopic interactions 

We consider granular assemblies made of nearly rigid 
equal-sized beads of diameter a and mass to, enclosed in 
a cuboidal simulation box. Beads interact through their 
contacts: the force transmitted is purely normal and is 
the sum of a Hertzian elastic term: 

= EV^h^/yS = ^KN{h)h, (1) 

and of a viscous term: 

F]:^ = C{mEy/^{ahy^^h ^ Cy/2mKN{h)h. (2) 

h is the normal elastic deflection, _B is a notation for 
E/{1 — i/^), where E is the Young modulus of the ma- 
terial the beads are made of, and v its Poisson ratio, 
and C is the level of viscous damping. K^ih) is the 
equivalent spring constant associated with the elastic 
force given by Eq. ([1]). Noteworthily, albeit nonlin- 
ear, Eqs H]) and ^ entail a velocity-independent nor- 
mal restitution coefficient eiqiC) in binary collisions. All 
the simulations reported here have been performed with 
C = 0.98 (cat = 3.3 X 10-3). This model has already 
been employed and discussed in several recent publica- 
tions [131, [l4| . Finally, as normal contact forces have no 
moment on spherical particles, their rotation is ignored. 



B. Boundary conditions and numerical tests 

Three different mechanical tests are numerically imple- 
mented to probe the solid behavior and the yield stress 
condition of the material. Those tests involve an external 
control on some of the entries of the Cauchy stress tensor 
CT. For a granular system at mechanical equilibrium, its 
expression involves the volume V of the system, the in- 
tergranular force Fij and the center-to-center vector 
for all pairs («, j) of contacting grains [2^[2^: 

^^■^^Fio® nj (3) 

Compressive stresses and shrinking strains are positive 
in our convention. 

In order to avoid any side wall effect, the simulation 
cell has periodic boundary conditions in all three di- 
rections (possibly affected by the Lees-Edwards proce- 
dure (30| when a non-diagonal stress component is im- 
posed). Simulation cell edges have lengths denoted as 
(-Z^a)i<Q<3 along the three coordinate directions of or- 
thonormal basis (e'Q.))i<Q.<3. Details on the equations 
governing the La 's and the possible shear strain variable 
may be found in Ref. • 

Before performing a mechanical test, an initial config- 
uration is prepared under isotropic pressure P with the 
same procedure as in jl^ [ij]. A granular gas of hard 
spheres, initially positioned on an FCC lattice, is ther- 
malized with collisions that preserve kinetic energy and 
then isotropically compressed [with the dissipative me- 
chanical model, Eqs ([I])-©] until a mechanical equilib- 
rium state is reached under presssure P. In the limit of 
small P, these isotropic equilibrated confi gur ations are 
the RCP states, as studied in Refs. [ni[ll,lli. 

Once prepared, the material may be subjected to var- 
ious loading paths. Three distinct quasistatic mechani- 
cal tests have been implemented, on externally applying 
stress tensor S: (i) Axisymmetric triaxial compression 
(TC) test: ^= Siei (g) ei + £26*2 (E) €2 + Use's <E) eg, 
with El = E2 < S3; (ii) Axisymmetric triaxial extension 
(TE) test : ^ = SiCi (g) ei + £262 (8)62 + Saea «) ea with 
El = E2 > S3; (in) Shear (S) test: ^ = P{ei (g) ei + 6*2 (?) 
62 -I- ea (g) ea) -I- T(ei (g) 62 -I- 62 (g) ei). 

Each test is employed to assess the material behavior in 
a particular direction of the principal stress space, which 
is the three-dimensional Euclidean space spanned by the 
stress tensor eigenvalues ci, CT2, and ca (the principal 
stresses). The principal stresses, if listed in decreasing 
order, are also denoted as tri > tru > am in the sequel. 
In equilibrium under the prescribed stress loading paths 
(TC, TE and S tests) their values are fisted in Tab. [T| 

In all implemented tests, pressure P = TrS/3 is kept 
constant while deviator stress S — PI is stepwise in- 
creased, with increments JS. We chose to apply iJSa — 
0.005 X P (or -0.005 x P) ifTTC (respectively: TE) tests, 
whence S^i = (5E2 = ±0.0025 x P, and St = 0.005 x P 
in S tests. In principal stress space, the load therefore 
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0-11 




TE 


(Til 


ail 


CT33 
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0"33 + |o"12| 


0"33 


o"33 — |cri2| 



initial state 



TABLE I: Principal stresses for the triaxial compression, tri- 
axial extension and shear tests. 




FIG. 1: (Color online) Sketch of the directions tested in a 
deviatoric plane with a triaxial compression test (with 0-3 = 
(Ti), a triaxial extension test (with 0-3 = om), and a shear test 
(with (cTi, 0-2,0-3) = (0-1,0111,011)). 



remains in a given deviatoric plane, i.e., a plane orthog- 
onal to the trisectrix ai =02 =0-3. The three studied 
stress paths are represented in Fig. [TJ 

For each prescribed stress tensor S, one -waits until a 
satisfactory mechanical equilibrium is reached before in- 
crementing S. A system is deemed equilibrated if the 
resultant force is zero on each bead, -with a tolerance 
set to lO^^a^P, and a^f^ — for each imposed stress 
component, with a relative error smaller than 10"'*. The 
calculation is stopped if the packing remains out of me- 
chanical equilibrium under the imposed stress tensor af- 
ter 5 X 10'' time steps and a total strain of 10%. The last 
value of S for -which an equilibrium state -was reached is 
kept as an estimate of the failure threshold. This proce- 
dure is schematized in Fig. [21 The same simulations are 
carried out on a number of different randomly assembled 
initial configurations to achieve statistical representativ- 
ity. Numerical results are averages over available sam- 
ples, the error bars sho-wn on the figures extending to 
one r.m.s. sample-to-sample deviation on each side of 
the mean value. 



C. Dimensionless control parameters 

Simulation results depend on a small set of dimension- 
less numbers, -which combine material properties and me- 




FIG. 2: Stepwise procedure employed to assess the failure 
properties of the material. 



chanical test parameters. Most definitions recalled belo-w 
are the same as in Refs. (Tsl, IT^. 
The stiffness parameter, defined as 




is such that n^^ measures the typical elastic deflection 
relative to particle diameter, h/a. Most simulations are 
conducted -with k, = 3.9 x 10'^ (corresponding to glass 
beads under P = 10 kPa [H [Tj). Half of the S tests 
have also been conducted -with k = 8.4 x 10'^. From p^ . 
-we kno-w that -with such stiffness levels, the difference 
bet-ween the various observables measured in the simula- 
tions and their values in the rigid limit of k — > +00 is 
smaller than the statistical uncertainty. 

Although some dissipation is necessary in the model to 
reach mechanical equilibrium, the level of viscous damp- 
ing C, is irrelevant in the quasistatic limit [l^ . Here 
it is set to 0.98, -whence a lo-w restitution coefficient, 
CAT = 3.3 X 10-3. 

When the cell is being deformed, -with strain rate e, 
the importance of inertia effects is characterized by the 
inertial number /, defined, as in Refs. [a,!!!, [13, ill, by 



The quasistatic limit corresponds to / — *■ 0. In order to 
avoid excessive acceleration of the system, a control on e 
is enforced, like in [13, US]) so that / never exceeds 10"*. 

Ratio 5E/P of deviator step to pressure is another con- 
trol parameter, -which should be kept to small values to 
track the quasistatic evolution of the system as accurately 
as possible. The values given in Sec. IHBI were observed 
to be satisfactory in this respect. As an example, TC 
tests with 5S3/P = 5 X 10"^ and ^Sa/P = 2 x 10"^ 
yield consistent results. 
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K I C 5T./P iV 

{8.4 X 10^3.9 X 10''} < 10"'' 0.98 0.005 {1372,4000,8788} 



TABLE II: Values taken by the dimensionless parameters. 

Finally, finite-size effects are expected [l3], hence a 
fifth dimensionless parameter in the problem, the number 
N of particles. Values of the dimensionless control pa- 
rameters in the presently reported simulations are listed 
in Tab. |lll We are chiefly interested in the macroscopic 
geometric limit, in which all mechanical properties are 
expected to depend on packing geometry alone, as an- 
nounced in the introduction. It was defined in |14l | as the 
triple limit of k ^ +oo (rigid limit), / ^ (quasistatic 
limit) and N —f +oo (thermodynamic limit). This limit 
was shown in Ref. [14| to be correctly approached with 
the range of parameters displayed in Tab. |TT1 
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TC 


4000 
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6.8° 


0.5° 


27.52 
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0.2° 


27.41 


0.03 




1372 
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8.6° 


0.5° 


27.81 


0.10 


TE 


4000 
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6.9° 


0.2° 


27.52 


0.02 




8788 
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6.0° 


0.4° 


27.40 


0.05 




1372 


6 


9.7° 


0.3° 


27.80 


0.04 
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4000 


10 


7.8° 


0.3° 


27.51 


0.07 




8788 


6 


7.0° 


0.4° 


27.41 


0.04 



TABLE III: Macroscopic friction angle (fi and Lade-Duncan 
parameter k, measured just before failure on 5'jv distinct ini- 
tial configurations, for different mechanical tests and different 
system sizes. Aip and Afe are the corresponding standard de- 
viations. 



III. FAILURE 

The material being initially assembled in an isotropic 
state, the range of stress tensors it will sustain in the 
solid state can be defined in principal stress space. From 
the known behavior of cohcsionlcss granular materials 
(with friction in the contacts) 0, [32, |33|, [S^ it is ex- 
pected - and it was explicitly checked in the case of S 
tests - that the boundary of the set of supported 
stresses is reached on increasing the deviatoric part of 
CT, away from the isotropic state. It is customary to de- 
fine a loading function (or yield function) / of principal 
stresses (cti, (T2, (T3), such that / < defines the region 
of supported stresses (which is believed to be convex in 
general [13]) and / = its boundary surface. 

For an assembly of perfectly rigid (k — +00), non- 
cohesive grains, the absence of stress scale implies that 
f{Xa_) = f{a_) for all A > 0. Thus, the failure surface 
of such a material has a conical shape in principal stress 
space. We assume this property to hold for our simulated 
system, which is close to the rigid limit. Consequently, 
it is sufficient to determine the intersection of the fail- 
ure surface with one deviatoric plane, that is the fail- 
ure curve. Furthermore, because of the isotropic prepa- 
ration method employed, / is a symmetric function of 
(ci, cr2, era) and the failure curve is left invariant by all 
permutations of (cti, (12, (T3). 

The (cohesionless) Mohr-Coulomb model 

/Mc(g) = o-i - (Tin - ((Ti + cTiii) sin if 

is often assumed (at least implicitly) true for granular 
materials ^35j. Numerous studies have been devoted to 
the macroscopic friction of sheared granular assemblies 
in various geometries [31[, and it is tempting to assume 
that the measured angle corresponds to friction angle (p 
in a Mohr-Coulomb model that would describe the failure 
properties of the material. The fourth column of Tab. lIIll 



displays the macroscopic friction angles measured with 
the three loading paths employed. It shows that <p de- 
pends on N, as already observed in Ref. but also on 
the kind of mechanical test employed. Consequently, the 
material cannot be described by a Mohr-Coulomb crite- 
rion. Although TC and TE tests are not sufficient to rule 
out the Mohr-Coulomb model since (p^^ — ip^'~^ is below 
the statistical uncertainties vitiating the results, the com- 
parison with the shear angles unambiguously invalidates 
this criterion. 

It would be appealing if one could characterize the fail- 
ure properties of the material with a single parameter 
that would not depend on the applied load direction. 
Such an attempt already proved successful for assem- 
blies of frictional equal-sized beads [3^, ji4] , whose failure 
curve was successfully modeled by a Lade-Duncan crite- 
rion ISal: 



jcri + o"II + Jiii)^ 
CTlCTiicriii 



(4) 



One should have fc > 27 in Q if condition / < is to 
define a non-empty cone of supported stresses (with k = 
27 only isotropic stresses would be possible in equilibrium 
and the material would behave like a liquid). According 
to Tab. mil this failure criterion also works well in the 
frictionless case: the values of k deduced from each of the 
three loading paths (see Tab. lIV[) agree with one another. 

Making use of the aforementioned permutation sym- 
metry, the stresses at failure computed for N = 1372 are 
plotted in Fig. [3l in the deviatoric plane, with the Lade- 
Duncan curve corresponding to fc(1372) and the three 
Mohr-Coulomb failure curves pertaining to the three dis- 
tinct numerical tests performed. The Lade-Duncan cri- 
terion is clearly the best model. 

Failure properties are, however, dependent on system 
size (Tab. HlH) . Fig. [5] plots the principal stresses at fail- 
ure in the deviatoric plane for TV = 1372, 4000, 8788 
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1 — sin ^ — sin^ + sin'^ ^ 
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1 + sin in — sin^ ip — sin^ ip 
27 



1 — sin ip 



TABLE IV; Mohr-Coulomb and Lade-Duncan parameters for 
the different loading paths, as functions of applied stress com- 
ponents, sintp, as defined in the second column, is used as an 
intermediate variable in the expression of parameter k in the 
third one. 




FIG. 3: (Color online) Calculated points with their error 
bars, Lade-Duncan criterion (k — 27.80, blue solid line) 
and cohesionless Mohr-Coulomb criteria corresponding to TC 
{ip = 8.3°, red dashed line), TE {p = 8.6°, red dotted line), 
and S tests {p> = 9.7°, red dotted and dashed line) for an 
assembly of A'' = 1372 particles. 




FIG. 4: (Color online) Principal stresses at failure and cor- 
responding Lade-Duncan fits for iV = 1372, 4000 and 8788 
(red dashed lines). The blue solid curve corresponds to the 
macroscopic limit (k — 27.22), whereas the two blue dotted 
curves, for k — 27.20 and k = 27.24, bound the uncertainty 
interval. 



A'^. Since other criteria predict a nearly circular failure 
curve for small deviatoric strength [ssj ] , other forms than 
the Lade-Duncan yield function could be fitted to the 
data in the macroscopic limit. First of all, we observed 
that the Driicker-Prager criterion [s^ (whose shape is al- 
ways circular in the deviatoric plane) does not correctly 
fit our results in the A^ +00 limit. The Matsuoka- 
Nakai criterion a model specifically tailored to cap- 
ture the failure properties of some sands, defined as 

J, , s (o"! + (Til + o-iii)((Ticrii -I- (Tiicriii -I- (Tiiicri) 
/mn(2:) = ^ m, 



(cri(TlI(Tin) 



with the corresponding Lade-Duncan fits. The domain 
bounded by the failure limit decreases with increasing A^. 
To evaluate the failure curve in the macroscopic limit 
of A'^ — *■ -|-c», principal stresses obtained in finite-size 
samples are extrapolated, assuming a linear dependence 
with A^~^/^. This assumption is proved to be statis- 
tically valid thanks to calculations [sj and the re- 
sulting principal stresses are plotted in Fig. 21 A Lade- 
Duncan fit of these extrapolated points with parameter 
= 27.22 ± 0.02 works well. As fcoo > 27, the failure 
surface is not reduced to the trisectrix in the N +00 
limit, and macroscopic systems can be equilibrated un- 
der moderately anisotropic loads, in agreement with p^ . 
The value taken by the Lade-Duncan parameter in the 
A' +00 limit corresponds to (^^^ = 4.4°±0.2° in triax- 
ial compression, (/?™ = 4.5° ± 0.3° in triaxial extension, 
and 1^%^ = 5.2° ± 0.3° for shear tests. The latter value 
agrees with the static friction angle given in Ref. [l4!| in 
the macroscopic geometric limit. 

One can observe that the shape of the criterion in the 
deviatoric plane becomes more rounded with increasing 



accurately fits the data extrapolated in the macroscopic 
limit (with m = 9.05 ± 0.01). Note however that this 
criterion is not suitable to describe failure in the smallest 
finite-size systems studied. 

In general, the expressions of failure criteria (Lade- 
Duncan, Matsuoka-Nakai, or Mohr-Coulomb) are purely 
phenomenological, and their justification is to provide 
a convenient fit function. In the present case, stress 
anisotropics will be related to other internal variables 
in Sec. IIVI but a prediction of the shape of the failure 
curve in the deviatoric plane (related to complex geo- 
metric properties of sphere packings) is currently beyond 
our reach. 



IV. SOLID BEHAVIOR AND 
MICROSTRUCTURE: THE ROAD TO FAILURE 

We now study the evolution of the material within the 
solid range, from the initial isotropic state to the fail- 
ure limit, with a particular emphasis on microscopic as- 
pects. We first investigate in Sec. IIV Al how the scalar 
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variables characterizing the internal state of the pack- 
ing evolve with growing deviator stress. Those vari- 
ables include solid fraction connectivity and coor- 
dination number, orientation-averaged pair correlations 
and force distributions, and were extensively studied 
in isotropic RCP states fl3 , IT3l |. Structural and force 
anisotropies [ill, ji^ . [H, l44j are studied in Sec. IIVBI In 
the spirit of [4^, we will show how stresses relate to 
anisotropy parameters. Then, Section IIV CI reports on 
the elastic moduli measured in nearly rigid anisotropic 
packings, with results generalizing previous numerical 
observations on isotropic RCP state elastic properties. 
Finally, the existence of a well-defined stress-strain law 
in the thermodynamic limit {N +00) is discussed in 
Sec. llVDl 



A. Scalar quantities 

The typical evolution of volume fraction $ with the 
deviatoric stress applied, characterized by sin if = (ui — 
o'ui) / {o'l + o'ui) , is depicted in Fig [51 It shows that what- 
ever the load applied, remains approximately equal 
to $RCP — 0.639 from the initial isotropic state to the 
failure threshold. In particular, the relative variations 
of $ remain smaller by more than an order of magni- 
tude than the deviatoric strains (see Sec. lIVPp . This 
is consistent with Ref. J^] which showed the material 
to be devoid of dilatancy in the macroscopic geometric 
limit. $ evolves quite erratically with simp; however, <i> 
seems to increase systematically when the applied stress 
is moderately anisotropic, then it reaches a maximum 
and finally, it decreases when the material approaches its 
failure limit. We have currently no convincing explana- 
tion for this phenomenon. The jumps in <i> are correlated 
to network rearrangements: we checked that the greater 
the jump, the more important the change in the contact 
list. 

The connectivity of the contact network is the set (p„) 
of probabilities for one grain to be involved in n con- 
tact forces. The coordination number z is linked to (p„) 
through z = "^ni^Pn- Average fractions Pn have been 
recorded for the three loading paths with N — 1372, 
4000 and 8788. At equilibrium, whatever the deviator 
applied, the set (p„) is found identical to the distribution 
measured on frictionless isotropic packings [T3j . For such 
packings, pi, p2 and ps vanish, because normal repulsive 
forces on a bead with less than four contacts cannot bal- 
ance. As in the case of isotropic packings, some grains, 
the rattlers, do not belong to the force-carrying structure: 
their proportion is estimated at po ~ 1.3%, which is close 
to the value obtained with isotropic packings 13]. In all 
simulations carried out with k — 3.9 x lO'*, the backbone 
coordination number z* — z{l — po)~^ remains equal to 
6.08±0.03 between the initial isotropic state and the fail- 
ure limit. By the isostaticity property of the backbone, 
z* tends toward 6 in the k +00 limit [Tl[. [T3. [l3] . 

If we now replace the contact network by network Ch 
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FIG. 5: (Color online) Volume fraction $ as a function of sin ip 
(as defined in Table for iV = 8788 and k = 3.9 x lO". Solid 
squares are for one TC test, solid triangles for one TE test 
and circles for one S test. Curves are stopped at the value of 
sin(/p corresponding to the failure limit. 



defined on declaring a bond to join all pairs of grains sep- 
arated by a distance smaller than h, then its coordination 
number z{h) is drawn as a function of h/a in Fig.[S]at the 
failure limit. Curves corresponding to the three studied 
loading paths are identical. z{h) starts from coordina- 
tion number z at h = and is the cumulated integral 
of the pair correlation function up to distance a + h be- 
tween sphere centers. One gets z{h) — z(0) oc {h/af'^ for 
h/a <^ 1 in all equilibrated packings. The same power 
law with exponent 0.6 has already been observed to fit 
z{h) data in the same range of gap h with isotropic pack- 
ings (RCP states) ^12, 13]. (No theoretical basis has been 
proposed for this power law, the prefactor and the expo- 
nent of which might slightly depend on the range of h 
fitted and on the treatment of rattlers [Tsl. [23t.) 

The probability distribution functions p{f) of normal- 
ized contact forces, / — F/jFY has a similar shape as 
reported in many numerical [H, [H, [H, [13, [11] and some 
experimental [i^, [s^] studies on granular media. p{f) 
first exhibits a slight increase, up to / ~ 0.5, and then 
it decreases, roughly exponentially for large /.Remark- 
ably, thanks to Kolmogorov-Smirnov tests [37[, we ob- 
served that all p.d.f. in the equilibrium configurations 
obtained for the different simulated stress states coin- 
cide. Neither the number of grains, nor the direction of 
the loading path, nor the proximity of the failure limit 
alter the force distributions, which remain statistically 
indistinguishable. p{f) thus coincides, within statistical 
uncertainties, with the form parametrized, e.g., in [l^ . 
As the backbone is isostatic, p{f) is geometrically deter- 
mined in the rigid limit. The p.d.f. may in particular be 
characterized by its moments, which we denote, for any 
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FIG. 6: (Color online) Average coordination number z{h) of 
network Ch as a function of h/a, computed from some equi- 
librated configurations near their failure limit under TC, TE 
and S tests. All three group of data collapse on a single curve. 
Inset: power law behavior of z{h) — z{0) for h/a <C 1, revealed 
by a double logarithmic plot. 



X > 0, as 



zix) = in 



in 



(5) 



and we obtain, e. g., Z{2) = 1.53 ± 0.02 and Z(5/3 ) 
1.29 ± 0.01 (those results will be useful in Sec. lT7U|) . 

Finally, using the same indicators as in [l3|, we ob- 
served no tendency towards the formation of locally crys- 
talline patterns in the configurations under varying devi- 
ator stresses. 



B. Anisotropy 

Previous works showed that the very origin of shear 
strength in granular materials is the anisotropy, both 
structural and mechanical, induced by the deviatoric 
stress dlj. We now explore this connection in the par- 
ticular case of frictionless, rigid bead assemblies. 

Mathematically, material anisotropy can be character- 
ized by the joint probability density function P(n, F) of 
finding an intergranular contact oriented along the unit 
vector n and carrying a force of intensity F. This quan- 
tity is of central importance since it intervenes in the 
expression of the Cauchy stress tensor. Bearing in mind 
that K ^ 1 and denoting the number of contacts by Nc, 



Eq. ([3]) can be rewritten as: 

Nr.a 



V 
V 



[F®n) 



j dn dFP{n, F)F' 
TV f 

= I (mE{n){F)fin®n 



niS) n 



(6) 



{F)f{ is the angular force density (it is equal to {F)/{Att) 
in the isotropic case) and E{n) is the probability density 
function of finding a contact along n. 



1. Structural anisotropy and its relation to stress ratios 

The anisotropy of the contact network is described by 
E{n). E is defined on the unit sphere of M-^, so it can 
be expanded in a series of spherical harmonics. Since 
contacts are undirected, odd order coefficients in the ex- 
pansion vanish. At the lowest order, the expansion is 
restricted to the spherical harmonics of order 2 and the 
coefficients are related to the second-order fabric tensor 
F = (n®fi). Furthermore, for shear tests, it was shown in 
Ref. |14| that a single anisotropic term of the expansion 
dominates: 



E{n) 



1 

47r 



Fi2 d^y{0,'ijj) 



(7) 



with dxy{9,Tp) = ISsin'^ 6'sin(2V;)/(87r) {9 is the colat- 
itude angle and tp the longitude angle of the spherical 
coordinates). In the case of a triaxial test, by axial sym- 
metry, the expansion of E in spherical harmonics up to 
the second order reads: 



33 



(8) 



with d^2 = 15(3 cos^ 6*- l)/(167r). 

Fig. [7] shows how the anisotropic term evolves with 
(CT33 — CTii)/(cr33-|-crii) (cth and (733 are principal stresses 
under a triaxial load), for systems of two different sizes 
subjected to a TC or TE test. Whatever the test per- 
formed, the absolute value of the anisotropic term in- 
creases with the applied deviator intensity. An analysis 
of the regression of fluctuations for the data of Fig. [7] in- 
dicates that the evolution of the anisotropic terms with 
stress deviator intensity tends to a well-defined curve in 
the macroscopic limit. The dependence is roughly lin- 
ear, even if one can notice that the slope of the curves 
seems to change around the boundaries of the solid range 
in the limit of iV ^ +00. (Other expressions involving 
principal stress ratios could have been used to charac- 
terize stress anisotropy). Although the maximum value 
of the anisotropy parameter is size-dependent, the slope 
Sfah of the straight line fitting in the macroscopic solid 
range the (sample-averaged) anisotropy parameter as a 
function of ((T33 — <Jii)/{<J^z + <Jii) for triaxial tests, and 
of (J\i/<Jti for shear tests, does not depend on N if the 
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FIG. 7: (Color online) Average evolution of the structural 
anisotropic term with r = (0-33 — an) /{ass + an) under TC 
(r > 0) and TE (r < 0) tests. Red crosses correspond to 
iV = 1372 and blue squares to iV = 8788. 



number of grains is large enough. For N > 4000, numer- 
ical simulation yield Sfab = 0.197 ± 0.010 for TC tests, 
S'fab = 0.210±0.015 for TE tests, and Sf^h = 0.158±0.015 
for S tests. 

The range of anisotropic pair correlations can be stud- 
ied by considering the fabric tensor of network Ch (defined 
in Section [IV Ap as a function of h. Anisotropy parame- 
ters are plotted as functions of h in Fig. [HI for maximum 
stress anisotropies (at the failure limit). They first de- 
crease for increasing h, and reach zero near h/a — 0.2. 
The small values of opposite sign measured at larger dis- 
tances are of the order of the statistical noise (~ 0.001) 
observed on isotropic configurations and should be in- 
terpreted with care. The spatial distribution of near, 
but distant neighbors thus tends to cancel the anisotropy 
of the distribution of contacting ones. The material 
anisotropy is short-ranged. In particular it is very nearly 
negligible on averaging over the complete first neighbor 
shell (i. e. up to the distance corresponding to the first 
minimum in t he p air correlation function, h/a ~ 0.35 
from Refs. 



2. Force anisotropy and its relation to stress ratios 

The mechanical anisotropy is described by the angular 
dependence of Like E(n), it can be expanded in 

a series of spherical harmonics of even order. To make 
things easier, only the expansion up to the second order 
is considered. As for the structural anisotropy, a single 
term, with the same symmetry, was assumed to domi- 
nate. For shear tests 
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FIG. 8: (Color online) Dominant structural anisotropic term 
at the failure limit as a function of the gap h ith = 8788 for 
TC tests (circles), TE tests (crosses) and S tests (triangles). 



and for triaxial tests 



{F)n - (^^ + ^33 d,- {6, V) j (F) (10) 

with (F) the average force intensity. 

Force anisotropy parameters H12 and H^^ are obtained 
by dividing the unit sphere in small regions. This allows 
to compute some values of (-F)a , and coefficients H12 and 
H23 are then derived by calculating the scalar product — 
defined as {f,g) = / (df^/(4^))/((?, V^), with / and 

g two functions defined on the unit sphere — of {F)fi with 
dxy and d^2 . 

The build-up of H^^, under a triaxial load for two differ- 
ent system sizes is displayed on Fig 121 It is very similar to 
the build-up of F32, — 1/3. The numerical data evidence a 
one-to-one correspondence with stress anisotropy, which 
is approximately linear for moderate deviators. The slope 
S'for of the plot of Fig. [S] seems to be independent of N 
when N is sufficiently large. With N > 4000, one has 
S'for = 0.250 ±0.012 for TC tests, S'for = 0.235 ±0.015 for 
TE tests, and Sfor = 0.173 ± 0.014 for S tests (for which 

H12 relates to stresses approximately as H12 = Sfor ). 

^22 



3. General connection between stress and anisotropy 

The observed relations between stress and fabric 
(Sec. lIVB ip or force (Sec. lIVB2p anisotropies were not, 
to our knowledge, previously reported in the literature. 
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FIG. 9; (Color online) Average evolution of the mechanical 
anisotropic term with r = ((T33 — crii)/((T33 + an) under TC 
(r > 0) and TE (r < 0) tests. Red crosses correspond to 
iV = 1372 and blue squares to iV = 8788. 



We argue below in Section IIVDI that they are specific to 
frictionless grains in the rigid hmit. 

Yet, a more general connection between stress and 
both fabric and force anisotropics can be derived on us- 
ing spherical harmonics expansions for the relevant stress 
components, as deduced from Eqs. (|8|10p for triaxial tests 
and from Eqs. (|7l9p for shear tests. Such a relation 
was repeatedly used for frictional systems, most often 
in 2D j4i,|4i,|5li. 

In the case of triaxial tests, keeping only the terms up 
to the second order yields 



E{n)x{F)n 



1 



-^33 + F: 



33 



1/3 



167r2 V An 
Combining this relation with Eqs 



0-11 



0'33 



V 

NMF) 
V 



127r Stt*-^^^ 



one gets 

- i^33 ' 1/3) 

- i^33 - 1/3) 



Consequently, one obtains 

"■33 ^^33 



33 



CTll 



(11) 

1 — H33 — F33 

In the case of shear tests, neglecting terms of order 
larger than 2 yields 



E{n) X {F)n 



1 



F12 + H12 



167r2 ■ V A-K 
By inserting the above equation in ([5]), one gets 
F12 + H12 NMF) 
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FIG. 10: (Color online) Numerical test of the approximations 
given by Eqs. dlTJ and ^ with iV = 8788 for TC tests 
(circles), TE tests (crosses) and S tests (triangles). 



hence the result: 



^ ^ HF,2 

022 



i?12). 



(12) 



Although Eqs. pT|) and (|12p are simple approximations, 
they work surprisingly well, as shown by Fig. 1101 

In Sections IIVB II and IIV B 21 fabric and force 
anisotropics were separately related to stress ratio (with 
approximate, linear relations involving parameters S'fab 
and iSfor)- Thus one should check for the consistency be- 
tween such observations and relations pT|) and In 
the case of triaxial tests, on writing down all quantities 
to first order in the small anisotropy parameters 7/33 and 
^33 — 1/3, one obtains the consistence condition: 

4 

•S'for + ftab = p:- (13) 

Similarly, for S tests one should have: 

^tor + ftab=i. (14) 

The values of S'fab and Sfor obtained in Sections IIV B II 
and IIV B 21 satisfy conditions (fT3|) and (I14p with good 
accuracy. 

The simple connection between stress and anisotropy 
parameters expressed by Eqs. (fTT|) and p2|l emphasizes 
the microscopic origin of a macroscopic quantity (a stress 
ratio in this case). In view of the different internal fab- 
ric symmetries in the triaxial and the shear tests, it is 
finally not surprising that the corresponding friction an- 
gles differ (whence the inadequacy of the Mohr-Coulomb 
criterion, Section ITTTl) . 
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In granular materials with friction, the shape of the 
particles influences the relative roles of geometry and me- 
chanics in the sustained stress [i^ . For frictionless spher- 
ical grains, near the failure limit, we find that the pa- 
rameters describing both anisotropics are approximately 
equal, so that half of stress ratio CT33/CT11 or cri2/o'22 is 
explained by geometric anisotropy and the other half by 
mechanical anisotropy. 

Despite those simple relations between stresses and 
anisotropy parameters, theoretically predicting the stress 
ratio at failure still remains a challenge. 



C. Elastic moduli 

The motivation for computing elastic moduli is 
twofold. First, elastic properties are usually more eas- 
ily measured in the laboratory than geometric data such 
as near neighbor correlations and coordination numbers, 
as discussed in Ref. . Then, the elastic moduli of fric- 
tionless bead packs under isotropic stresses were stud- 
ied by numerical simulations [ll|, [l^l , and shown to ex- 
hibit singular properties, which we now seek to generalize 
to anisotropic stress states. Specifically, while the bulk 
modulus, B, shows little difference with well coordinated 
frictional packings [g^l, the shear modulus, G, is anoma- 
lously small. G/B tends to vary proportionally to the 
degree of force indeterminacy ^24, 52], which vanishes in 
the rigid limit, as k^^^^. Isotropic frictionless bead packs 
also possess stiffness matrices (or "dynamical matrices" ) 
with an anomalous distribution of eigenmode frequen- 
cies 11], which stems from the nearly isostatic character 
of the contact network [ssj . 

For simplicity, we restrict our investigations to the elas- 
tic moduli of equilibrium configurations obtained in TC 
or TE tests. They are numerically evaluated on building 
the stiffness matrix of contact networks and solving lin- 
ear systems of equations for displacements in response to 
small load increments, as explained in The results 
are devoid of size effects and sample to sample fluctua- 
tions regress as N increases. There are five independent 
elastic constants in such cases (a number which would in- 
crease to nine for simple shear tests), which express a lin- 
ear relation between stress increments A.aij and strains 
Eij, from a reference equilibrium anisotropic state, as 
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(15) 

The material symmetries - invariance by rotation around 
axis 3 and by symmetry about all three planes of coordi- 
nates - determine the form of the matrix of elastic moduli 
in (fT5|) . and also request that Cu — C12 = 2C55 (express- 



ing the equality of two shear moduli in plane 1,2). Such 
symmetries are very well satisfied (one has, e.g., C13 = 
C23 with relative errors smaller than lO^'^ for N — 8788). 
The moduli in the initial isotropic state all relate to B 
and G as Cn = S + 4G/3 = C33, G12 = B-2G/3 = C13, 
G44 = G55 = G. Then, longitudinal moduli (i.e., Cu, 
with 2, 3) are larger in the direction of the major 
principal stress: thus one observes C33 > Cn in triaxial 
compression and the opposite inequality in triaxial ex- 
tension. This corresponds to different longitudinal sound 
wave velocities Cu/ pm {pm denoting the mass density 
of the material) propagating in direction 3 and in the or- 
thogonal plane. Such anisotropies of the elastic moduli 
were reported in the literature on sands [H, [H, [sH] and 
bead packings fs?]. They can be attributed to the ef- 
fect of both anisotropies, of fabric and forces, evidenced 
in Section IIVBI the material is stiffer in the principal 
stress direction because it is favored in the distribution 
of contact orientations, and also because contacts nearly 
parallel to this direction tend to carry larger forces. As 

1 /3 

Hertz's law, Eq. ([1]), entails that Kj^ cx FjJ , such con- 
tacts are stiffer. To sort out the possible effects of fabric 
and force anisotropies, we computed elastic moduli both 
for the Hertzian contact model and for linear contact 
elasticity, with some constant, force-independent contact 
stiffness Kn- (As the packing geometry is very nearly 
that of a set of rigid beads, statistically similar config- 
urations would have been obtained on simulating bead 
assemblies with linear unilateral elastic contact forces). 
We focus in the sequel on the upper left square block of 
order 3 within the matrix of moduli written in Eq. (llSp . 
which we denote as c. All of its elements are larger by 
about 2 orders of magnitude than shear moduli C44 and 
G55 , whatever the stress anisotropy. The ratio of all other 
elements of matrix c to C33 are plotted in Fig. [Tl] for 
Hertzian and for linear contact elasticity. The variations 
of C11/C33 with cr33/crii shown on Fig. [TT] are qualita- 
tively expected. More surprisingly, since the moduli eval- 
uated with linear contact elasticity are not sensitive to 
force anisotropy, the dependence of such ratios on stress 
anisotropy is about the same for both contact laws. 

Such results, as we now explain, are due to the pecu- 
liar nature of matrix c. Let si, §2, S3 denote unit vectors 
in the space of stress or strain tensors with eigendirec- 
tions parallel to the coordinate directions, forming an or- 
thonormal basis in which coordinates are Acru (i=l, 2, 3) 
(or eu) for stress (resp. strain) increments. Matrix c de- 
fines a linear operator within this space. Under isotropic 
stresses, c has eigenvalues Cj = 3B, Cu — Cm = 2G, 

and eigenvectors are §1 = (si + §2 + S3)/V3, and any 
pair of vectors orthogonal to Si . The increment of stress 
in direction Si is proportional to the preexisting equilib- 
rium stress tensor (denoted here as vector P\/3Si). As 
an approximation, since G/ 3> G// and G/ ^ Cm, one 
may write: 

£~ G/S"! (gjS*!, (16) 

bearing in mind that the right-hand-side is of course a 
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FIG. 11: (Color online) Ratios C11/C33 (red), C13/C33 
(pink), and C12/C33 (blue), in TE and TC tests, versus princi- 
pal stress ratio, for Hertzian (square dots) and linear (round 
dots) contact elasticity and A'^ — 8788. Continuous (black) 
lines correspond to the predictions of Eqs. (|18|) . 



FIG. 12: (Color online) Dominant eigenvalue Ci of tensor of 
elastic moduli, for Hertzian and linear contact elasticity (A'' — 
8788), versus principal stress ratio along TE and TC loading 
paths, compared to the predictions of Eqs. (|19|l and ((20}, 
depicted, due to the slight statistical uncertainty, as narrow 
zones between horizontal dashed lines. 



singular matrix. On using (|16[) . all moduli would be 
equal to B in the isotropic state and all ratios equal to 1 
in Fig.[Tl]for an = (733, which is very nearly satisfied. In 
[24] , it was argued that the "dominant" modulus, B is in- 
sensitive to the "barely rigid" character of the nearly iso- 
static contact network because it expresses the response 
to a load increment proportional to the preexisting load. 
We now apply similar ideas to anisotropic stress states. 
We first define Si as the unit vector proportional to the 
preexisting, equilibrium stress. The loading parameter 
in triaxial loading paths may be defined as a such that 
CTii = 0-22 = (1 — a)P while 0-33 = (1 -I- 2a)P. We thus 
set: 



Si = 



1 



V3 + 6a' 



[(l-a)(si+S2) + (l + 2a)s3]. (17) 



We observed 6*1 to be, with very good approximation, an 
eigenvector of c, with eigenvalue Ci close to its value in 
the isotropic state. Due to the material symmetries in TC 
and TE tests, the second eigenvector should be S2 = {§1 — 
S2)/V^, a property also well satisfied by the numerical 
data - and the third one is of course orthogonal to 5*1 and 
52- We observed the corresponding eigenvalues C// and 
Cm to remain below 0.02 x C/ in all cases, whatever 
the stress anisotropy and the contact law (Hertzian or 
linear). Thus it is possible to approximate matrix c on 
using relation p^ . with definition (fTT]) for vector Si. 
This yields theoretical expressions for the ratios between 
moduli: 



C33 



C12 

C33 ' 



C33 



1 



l + 2a 



(18) 



Fig. [TT] shows that those approximations are quite ac- 
curate. Thus stress anisotropics influence the tensor of 



elastic moduli in a peculiar way, due to its nearly uniax- 
ial, singular structure, which is independent of the con- 
tact law. In a good approximation all moduli, except the 
very small, singular ones, are proportional to C/ with 
coefficients that are determined by the stress state. 

On exploiting the isostaticity property of the contact 
network, it turns out that the dominant eigenvalue of 
tensor c, C/ can be written, in very good accuracy, as 
a simple function of solid fraction $, coordination num- 
ber z and moments of the (geometrically determined) 
force distribution. Such a relation was established for 
the bulk modulus B of isotropic states in [13] , where it is 
called the Reuss estimate. In general, it provides a lower 
bound to the modulus, which becomes exact when force 
increments are proportional to preexisting forces. This 
condition is exactly fulfilled by the response of isostatic 
contact networks to an increment of stress tensor that is 
proportional to the preexisting stress tensor. On adapt- 
ing the approach followed in [24] to the case of anisotropic 
stress states, one readily obtains, in the case of Hertzian 
contacts: 



Ci = Cf = 



31/3 



2Z(5/3) 



TT 



2/3 



(19) 



For linear contact elasticity, the corresponding prediction 
reads 



L _ z^Kn 



(20) 



Z(5/3) and Z{2) values are given after Eq. [5J All quan- 
tities appearing in those formulas were observed in Sec- 
tion IIV Al to remain constant throughout the range of 
supported stresses. Thus C/ should not depend on prin- 
cipal stress ratio. Fig. [T^ shows that the numerical data 
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abide very well by the predictions of Eqs ((T9)) and ([20)l . 
Thus, all moduli, except the soft ones that vanish in the 
rigid limit, are predicted. In the case of simple shear, we 
expect similar properties to apply, on adequately redefin- 
ing 5*1 in the direction of the applied load. In general, 
for arbitrary applied stresses within the supported range 
/(ct) < defined in Section HITl one should have a nearly 
uniaxial tensor of elastic moduli. 

The stress increment or strain range for elastic re- 
sponse is expected to shrink to naught in the double limit 
of K — > -|-oo and iV — s- cx), like the stability range of a con- 
tact network 20]. Thus, in practice, in order to observe 
the peculiar elastic properties of nearly rigid frictionless 
bead assemblies, one should adequately choose stiffness 
level K, which should be large enough to approach the 
rigid limit but small enough for some elastic response to 
be measurable. Interestingly, poorly coordinated pack- 
ings of frictional disks [ssl, or spheres [U, [gOl tend to 
exhibit similar elastic anomalies, although, most often, 
in a weakened form, because such systems do not spon- 
taneously form isostatic contact structures in the rigid 
limit [13| . Even though truly frictionless particles do not 
exist in the laboratory, our results might therefore bear 
some relevance in more general situations of contact net- 
works with quite a small level of force indeterminacy. 



D. Constitutive relations 




FIG. 13: (Color online) Evolution of strain £33 with ass/cru 
in triaxial tests, for A'^ — 1372 (red crosses connected by a dot- 
ted line), A'^ — 4000 (brown triangles connected by a dashed 
line), and TV = 8788 (blue squares connected by a solid line). 
Results are averaged over all available samples, and restricted 
to the macroscopic solid range. 



The observations of Sections IIVBI and IIV CI are 
strongly reminiscent of the results obtained on dealing 
with exactly rigid frictionless grains [l^, [6l| . If contacts 
are rigid, the response to an applied stress increment pro- 
portional to the preexisting stress is also perfectly rigid: 
the corresponding strain is exactly zero. For all other 
stress increment orientations, the rigid contact network, 
in the limit of iV 00, has to rearrange The 
resulting strain is determined by the geometry of the 
packing, rather than by some material stiffness. Thus, 
using the notations of Sec. IIV C[ C/ is infinite, while 
Cii = Cm = 0. This behavior also entails a one-to-one 
correspondence between stress and fabric anisotropy, in 
agreement with Sec. IIV B II Because of isostaticity, the 
force distribution is completely determined by the force 
network, whence the relation evidenced in Sec. IIVB 21 
In this respect, assemblies of frictionless grains differ 
from systems with intergranular friction, in which one 
given contact network may support stresses within a fi- 
nite range in the thermodynamic limit, for arbitrary large 
stiffness levels k [61| (whence vertical parts in stress ver- 
sus strain plots, as obtained in simulations with models 
of rigid grains [5l|, [H, [G^I). This property of frictional 
grain assemblies excludes the possibility of a one-to-one 
relation between stresses and fabric. 

Rigid frictionless grain assemblies, on the other hand, 
were reported [l^, [6l| to be devoid of the stress-strain 
relations (which depend on loading history) obtained in 
simulations of model frictional systems [33, [H, HH , and 



classically modeled, for sands, in soil mechanics [l^. [l7l|. 
This conclusion was based on a statistical analysis of the 
strain response to stress increments, which was modeled 
as a Levy-distributed random variable 6J|, precluding 
the regression of strain fluctuations in the thermody- 
namic limit. Such results contrast with the ones obtained 
with particles interacting with soft potentials, such as 
Lennard-Jones glasses, in which case fluctuations around 
the average stress-strain curve were explicitly shown to 
regress in the thermodynamic limit (65} . 

In the present case, the macroscopic mechanical re- 
sponse is also dominated by packing rearrangements: 
macroscopic strains are much larger than typical contact 
deflections (of order k~^). Macroscopic strains, as plot- 
ted versus applied stress ratio along the triaxial test paths 
in Fig.[T21 do not appear to behave like a Levy flight tra- 
jectory: results pertaining to the two larger sample sizes 
tend to cluster around the same average curve. However, 
the regression of fluctuations in the limit of A/ 00 is 
much less clearcut than in the results of, e.g., Fig. [71 er- 
ror bars are only very slightly reduced between N = 1372 
and N — 8788, and still extend to a notable fraction of 
averages (typically 30%). Our data very likely provide 
unsufRcient statistics because of sample size limitations, 
and larger systems should be studied. Yet it is tempting 
to speculate that large enough samples, for given k, do 
approach a well-deflned stress-strain behavior for given 
loading paths, but that their size should exceed a cer- 
tain characteristic length ^ that diverges in the limit of 
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K — > +00. In this interpretation, for any given value of 
K, samples of (linear) size below £ would exhibit the sin- 
gular behavior observed in Ref. [20| (in which rigid con- 
tacts were simulated, with a specific numerical technique 
exploiting the isostaticity property). Only for samples 
larger than ^ (and hence, for larger and larger samples 
as K is increased) should one recover a well-defined stress- 
strain relationship in monotonic loading. Further inves- 
tigations of this conjecture are beyond the scope of the 
present paper. 

V. DISCUSSION 

The present study generalizes the results on the macro- 
scopic friction of frictionless bead packs, previously ob- 
tained in simple shear, to other loading paths, and pro- 
poses a form of the failure criterion valid for arbitrary 
stress directions. This failure condition is somewhat dif- 
ferent from the Mohr-Coulomb condition and best ex- 
pressed in the Lade-Duncan form. As previously ob- 
served despite rather strong finite size effects, the 
system is able to sustain finite stress deviators in the 
macroscopic geometric limit, in which the Lade-Duncan 
parameter, evaluated at koo — 27.22 ± 0.02, is to be re- 
garded as a basic geometric property of disordered sphere 
assemblies. Changes of volume fraction $ as deviatoric 
stresses evolve from zero to yield threshold values are 
quite small and erratic (in spite of a very slight tendency 
toward contractance under small deviator, and to volume 
increase close to failure) and might be neglected, given 
statistical uncertainties, in a first approach. Thus $ re- 
mains approximately equal to the RCP value. All classi- 
cal characterizations of packing geometry and force net- 
works by scalar or orientation- averaged variables, includ- 
ing the distribution of normal forces, do not distinguish 
anisotropic equilibrium states from the initial isotropic 
structures equilibrated under hydrostatic pressure. 

Thus the equilibrated configurations may be regarded 
as anisotropic random close-packing states. Isotropic 
RCP states, in the limit of rigid particles, are local min- 
ima of sample volume in configuration space, under the 
constraint of impenetrability of particles. Anisotropic 
ones also minimize the potential energy of the applied 
stresses, viz. 

where strain tensor e, assumed small, has to be defined 
with respect to some arbitrary reference configuration. 
Consequently, they do not maximize volume fraction $, 
and, although stable equilibrium states, do not qual- 
ify as "strictly jammed" according to the definition of 
Refs. [6^ [131 ■ That their volume fraction is no smaller 
(and occasionally slightly larger) than $rcp obtained in 
isotropic configurations is due to the multiplicity of differ- 
ent possible equilibrium networks and minima of poten- 



tial energies W, which are not connected by quasistatic 
trajectories. 

Fabric and force anisotropics can be efficiently charac- 
terized with one coefficient in an expansion in spherical 
harmonics. Each one of such coefficients is a function of 
stress anisotropy. The existence of such relations is spe- 
cific to frictionless systems, in which any change of stress 
direction tends to entail rearrangements and changes in 
the contact network. Meanwhile, like in granular systems 
with friction, stresses can be expressed, in good approx- 
imation, as combinations of fabric and force anisotropy 
parameters. 

Elastic moduli exhibit similar anomalies in the rigid 
limit as in isotropic states, with a nearly uniaxial tensor 
of elastic moduli, the dominant eigenvalue of which (the 
only non-singular one) expresses the response to load in- 
crements parallel to the preexisting load in stress space. 
Meanwhile, the moduli in orthogonal directions vanish in 
the rigid limit, as in isotropic systems (and the "density 
of states" for eigenmodes is expected to exhibit the same 
singularities |53[). These properties can be expected to 
apply to any situation of very small force indeterminacy 
in particle packings. 

Our results seem to indicate that a deterministic stress- 
strain curve for monotonic loading along given devia- 
toric paths should be obtained in the macroscopic limit, 
thereby contradicting the conclusions of [lO] , based on an 
exactly rigid system in 2D, although the simulated sam- 
ples still seem too small to reach a clear conclusion about 
the regression of strain fluctuations for given applied 
stresses. This point obviously deserves further investiga- 
tions, as well as the spatial structure and displacement 
correlations in deformation and rearrangement mecha- 
nisms. The possibility of a diverging length scale in the 
rigid limit of At — > -l-oo (entailing the non-commutation 
of the limits of k — > -l-oo and of iV — > +00) should be 
explored in further simulations of larger systems with 
varying stiffness level. 

Another issue worth investigating is that of the pos- 
sible uniqueness of equilibrium states, in the statistical 
sense, under a given supported state of stress. Just like 
simulation results appear to support the idea of a unique 
RCP state under isotropic pressure [l3|, provided a fast 
enough assembling process bypasses crystal nucleation, 
the results reported here suggest that the internal state of 
the packing in equilibrium could be uniquely determined 
by the current value of stresses, whatever the loading 
history. Such a conjecture is, in particular, supported by 
the observation of a one-to-one correspondence between 
stress and all measured internal state variables, such as 
fabric or force anisotropy parameters. 

Eventually, we expect that the knowledge of the be- 
havior of frictionless granular assemblies will be useful in 
the design of compaction strategies (lubrication, vibra- 
tion, cyclic loading...), which can be regarded as meth- 
ods to circumvent the influence of friction [l^. Other 
interesting perspectives involve the treatment of differ- 
ent particle shapes [H, and polydispersities [i^ . 
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